/*
 *  pulse.cpp laser pulses train differents characteristics by pulse
 *
 *  Created by Alexis Chacón and Camilo Ruiz
 *I0
 */

#ifndef LASER_H
#define LASER_H
#include <stdlib.h>
#include "timegrid1.h"
#include "timeobject1.h"
#include "constant.h"
#include<vector>
using namespace std;

class laser 
{

public:

	int npulses;
	int nmaxt;

	double tstart_first;
	double dt;
	double offset;
	double outset;		

	double major0;
	double minus0;

	vector<double> tstart;
	vector<double> thalf;
	vector<double> tend;

	vector<double> I0;
	vector<double> e;	
	vector<double> E0x;
	vector<double> E0y;
	
	vector<double> w0;
	vector<double> period0;
	vector<double> cycles0;
	
	vector<double> cep0;
	vector<double> phi_rel;	
	vector<double> delay0;
	
	timegrid1 g;
	timeobject1 efield;
	timeobject1 avector;

	timeobject1 *ef;
	timeobject1 *env;
	timeobject1 *av;

	timeobject1 *aux1;
	timeobject1 *aux2;
	
	/*==========================*/
	/*      MAIN FUNCTIONS
	 /*==========================*/
	laser(int _npulses);      							                                //Creator Object
	void laser_pulses(double _dt, double _tstart_first, double _offset, double _outset);	//LASER PULSES TRAIN
	
	/*=========================*/
	/*   SECUNDARY FUNCTIONS
	 /*=========================*/
	inline void Initialize_AmplitudeAndPeriod();       //Function initialize max amplitude and period
	inline void Set_StatAndEndTime();		           //Function "Start" and "end" time by each pulse
	void Laser_Grid();		                           //Function object time axis and field_set
	void Evaluation_Laser_Pulse();	                   //Funtion evaluation pulses train
	void Sum_Pulses();		                           //Function pulses sum
	void set_VectorPot();		                       //Function vector potential by each pulse
	void vec_potential();	                           //Function vector Potential to pulses train

	double qmajor(vector<double>& v);
	double qminus(vector<double>& v);
};


//==============================================================================//
			/*=== MAIN FUNCTIONS ===*/

/*===========================================================          		
		/*=== OBJECT'S CONSTRUCTOR LASER  ===*/
laser::laser(int _npulses)
  {
	  npulses=_npulses;
	  tstart.resize(npulses,0.0);
	  thalf.resize(npulses,0.0);	  
	  tend.resize(npulses,0.0);   
     
	  I0.resize(npulses,0.0);
	  e.resize(npulses,0.0);	  

	  E0x.resize(npulses,0.0);
	  E0y.resize(npulses,0.0);	  
     
	  w0.resize(npulses,0.0);
	  period0.resize(npulses,0.0);
     
	  cycles0.resize(npulses,0.0); 
	  cep0.resize(npulses,0.0);
	  
	  phi_rel.resize(npulses,0.0);
	  
	  ef  = new timeobject1[npulses]; //reserve memory
	  env = new timeobject1[npulses]; //reserve memory
	  av  = new timeobject1[npulses];
	  aux1  = new timeobject1[npulses];	  
	  aux2  = new timeobject1[npulses];	  
     
	  if (npulses > 1 )
		  delay0.resize(npulses-1,0.0);	  		  
  }//End initialize variable  


/*===========================================================          		
 /*===  LASER PULSES FUNCTION  ===*/
void laser::laser_pulses(double _dt, double _tstart_first, double _offset, double _outset)
{
   	/*======= Pulse's Parameter ======*/
   	tstart_first=_tstart_first;                     //Start time first pulse
   	dt=abs(_dt);			                        //Time increment 	
   	offset=abs(_offset);	                        //Offset left of the pulse train
   	outset=abs(_outset);	                        //Offset left of the pulse train	
	
   	//================================//
	Initialize_AmplitudeAndPeriod(); 		     	//Initialize max amplitude and period        
   	Set_StatAndEndTime();	 	     	            //"Start" and "end" time by each pulse
	
   	major0=qmajor(tend);	     	                //Maximum time to all pulse
   	minus0=qminus(tstart);          	            //Minimun time to all pulse     
	
   	Laser_Grid();       	     	                //Objet time axis      
   	Evaluation_Laser_Pulse();          	            //Evaluation pulses
	
   	Sum_Pulses();	     	     	                //Pulses sum (build pulses train)
   	set_VectorPot();		     	                //Vector potential by each pulse
   	vec_potential();	     	                    //Vector potential to pulses train 

}//End laser_pulses

/*============ END MAIN FUNCTIONS ===============*/



/*=================================================================*/
		/*========= SECUNDARY FUNCTIONS ===========*/

//== Function initialize max amplitude and period ==// 
inline void laser::Initialize_AmplitudeAndPeriod()
{
	
	for (int i=0;i<npulses;i++)
	{
		double factor = 1.0/(1.0+e[i]*e[i]);          //The ellipticity is e = E0x/E0y
		E0x[i]        = e[i]*sqrt(I0[i]*factor/3.5e16);
		E0y[i]        = sqrt(I0[i]*factor/3.5e16); 		
		period0[i]    = dospi/w0[i];
	}
}


//== Function "start" and "end" time by each pulse ==//
inline void laser::Set_StatAndEndTime()
{
	for (int kpulse=0;kpulse<npulses;kpulse++)
	{		
		if(kpulse==0)
		{
			tstart[kpulse] = tstart_first;
			thalf[kpulse]  = tstart[kpulse] +cycles0[kpulse]*period0[kpulse]/2.0;
			tend[kpulse]   = tstart[kpulse]+cycles0[kpulse]*period0[kpulse];
		}
		else
		{
			tstart[kpulse] = tend[kpulse-1] + delay0[kpulse-1]-(cycles0[kpulse-1]*period0[kpulse-1] + cycles0[kpulse]*period0[kpulse])/2.0;
			thalf[kpulse]  = tstart[kpulse] + cycles0[kpulse]*period0[kpulse]/2.0;
			tend[kpulse]   = tstart[kpulse] + cycles0[kpulse]*period0[kpulse];
		} 
	}//End loop
}


//== Function object time axis and field_set  ==//
void laser::Laser_Grid()   
{	
	nmaxt=floor((major0+outset-(minus0-offset) )/dt)+1;
	g.set_grid(nmaxt, dt, minus0-offset);      
	efield.put_on_grid(g);
	avector.put_on_grid(g);
	
	for (int kfield=0; kfield <npulses;kfield++){
		env[kfield].put_on_grid(g);
		ef[kfield].put_on_grid(g);
		av[kfield].put_on_grid(g);
		aux1[kfield].put_on_grid(g);
		aux2[kfield].put_on_grid(g);		
	}
	cout << "Ntime: "<<nmaxt<<endl;
	cout << "Memory: "<<((2+3*npulses)*24+8)*nmaxt*1e-6 << "  Mb"<<endl;
}


//== Funtion evaluation pulses train ==//
void laser::Evaluation_Laser_Pulse()
{
	for(int kpulse=0;kpulse<npulses;kpulse++)
	{
		for(int ktime=0;ktime<nmaxt;ktime++)
		{
			double arg0=w0[kpulse]*(g.t[ktime]-tstart[kpulse])/2/cycles0[kpulse];            //Envelope Argument
			
			double arg1=w0[kpulse]*(g.t[ktime]-tstart[kpulse])+cep0[kpulse];                 //Carrier Phase x
			double arg2=w0[kpulse]*(g.t[ktime]-tstart[kpulse])+cep0[kpulse]+phi_rel[kpulse]; //Carrier Phase y
			
			if (g.t[ktime] < tstart[kpulse] || g.t[ktime] > tend[kpulse]){
				ef[kpulse].f[ktime][0]  = 0.0;
				ef[kpulse].f[ktime][1]  = 0.0;
				
				env[kpulse].f[ktime][0] = 0.0;
				av[kpulse].f[ktime][1]  = 0.0;
			} 
			else{
				env[kpulse].f[ktime][0] = E0x[kpulse]*sin(arg0)*sin(arg0);			   
				env[kpulse].f[ktime][1] = E0y[kpulse]*sin(arg0)*sin(arg0);			   
				
				ef[kpulse].f[ktime][0] = env[kpulse].f[ktime][0]*sin(arg1);			   
				ef[kpulse].f[ktime][1] = env[kpulse].f[ktime][1]*sin(arg2);
			}
			av[kpulse].f[ktime][0]=ef[kpulse].f[ktime][0];
			av[kpulse].f[ktime][1]=ef[kpulse].f[ktime][1];		   
	   }
	}
}


//== Function pulses sum (build pulses train) ==//
void laser::Sum_Pulses()
{
	//Start loop sum pulses
	for(int ktime=0;ktime<nmaxt;ktime++)
	{
		efield.f[ktime][0]=0.0;
		efield.f[ktime][1]=0.0;		
		for(int kpulse=0;kpulse<npulses;kpulse++)
		{
			efield.f[ktime][0]+= ef[kpulse].f[ktime][0];
			efield.f[ktime][1]+= ef[kpulse].f[ktime][1];

		}	
	}//End sum set pulses       
}



//== Function vector potential by each pulse ==//
void laser::set_VectorPot()
{
	for (int kpulse=0;kpulse<npulses;kpulse++)	
		for (int ktime=0; ktime<nmaxt; ktime++){
			aux1[kpulse].f[ktime][0]=av[kpulse].f[ktime][0];
			aux2[kpulse].f[ktime][0]=av[kpulse].f[ktime][1];					
		}
	
	for(int kpulse=0;kpulse<npulses;kpulse++)
	{
		aux1[kpulse].integrate();
		aux2[kpulse].integrate();		
		double aa=aux1[kpulse].f[0][0];
		double bb=aux2[kpulse].f[0][0];
		
		for(int ktime=0;ktime<nmaxt;ktime++)
		{
			av[kpulse].f[ktime][0]=aux1[kpulse].f[ktime][0]-aa; 
			av[kpulse].f[ktime][1]=aux2[kpulse].f[ktime][0]-bb;
			
			av[kpulse].f[ktime][0]*=-1; 	//According to Landau & Lifshitz in The Classical Theory of Fields,
			av[kpulse].f[ktime][1]*=-1; 	//and Remetter Attosecond Wave Packet Interference, vector potential has negative sign.
			
			/*if (g.t[ktime] < tstart[kpulse] || g.t[ktime] > tend[kpulse]){
				av[kpulse].f[ktime][0]=0.0; 
				av[kpulse].f[ktime][1]=0.0;
			}*/
		}  
	}
}


//== Function vector potential to pulses train ==//
void laser::vec_potential()
{
	for(int ktime=0;ktime<nmaxt;ktime++)
	{
		avector.f[ktime][0]=0.0;
		avector.f[ktime][1]=0.0;		
		for(int kpulse=0;kpulse<npulses;kpulse++)
		{
			avector.f[ktime][0]+= av[kpulse].f[ktime][0];
			avector.f[ktime][1]+= av[kpulse].f[ktime][1];
		}	
	}//End sum set pulses       
}


//Start function mayor
double laser::qmajor(vector<double>& v)
{
	int n=v.size();
	double may = v[0];
	
	for (int k=1;k<n;k++)   
		if (v[k]>may) may=v[k];
	
	return may;
}//End mayor

//Start function minus
double laser::qminus(vector<double>& v)
{
	int n=v.size();
	double min=v[0];
	
	for (int k=1;k<n;k++)
		if (v[k]<min) min=v[k];	
	
	return min;
}//End minus

/*========================== END  SECUNDARY FUNCTIONS =============================*/

#endif
//END
